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Abstract. The SD-SPIDER method for the characterization of ultrashort laser 
pulses requires the solution of a nonlinear integral equation of autoconvolution 
type with a device-based kernel function. Taking into account the analytical 
background of a variational regularization approach for solving the correspond¬ 
ing ill-posed operator equation formulated in complex-valued L^-spaces over fi¬ 
nite real intervals, we suggest and evaluate numerical procedures using NURBS 
and the TIGRA method for calculating the regularized solutions in a stable 
manner. In this context, besides the complex deautoconvolution problem with 
noisy but full data, a phase retrieval problem is introduced which adapts to the 
experimental state of the art in laser optics. For the treatment of this problem 
facet, which is formulated as a tensor product operator equation, we derive well- 
posedness of variational regularization methods. Case studies with synthetic and 
real optical data show the capability of the implemented approach as well as its 
limitation due to measurement deficits. 
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1. Introduction 


About two decades ago and motivated by a problem from spectroscopy, namely 
the evaluation of highly resolved functions of the density of unoccupied states 
from appearance-potential spectra (cf, e.g., [291 ), the inverse problem of deauto¬ 
convolution came into the focus of the mathematical literature for the first time 
(cf. [21116]). In particular, the stable approximate determination of real functions 
/ over M with compact support, say supp(/) C [0,1], from noisy data of its self¬ 
convolution g with suppl^f) C [0, 2] became of interest, which is equivalent to the 
solution of the Volterra type nonlinear integral equation 


S 


( 1 ) 



0 < s < 2, 


0 
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occurring also in stochastics if the probability density function /(r) with support in 
[0,1] of a random variable j£ is to be determined from data of the density function 
g{s). In this context, g{s) corresponds to the random variable 2) + 3, where 
and 3 are stochastically independent random variables and 1^,3 are identically 
distributed. It was shown in [TB] resp. HH that Q written as an operator equation 
with the nonlinear forward operator of autoconvolution between the real spaces 
L^(0,1) and h^(0,1) resp. L^(0,2) of quadratically integrable functions is locally 
ill-posed everywhere in the sense of m Dehnition 2]. Furthermore, regularization 
approaches were developed and evaluated, where the astonishing fact appeared 
that in spite of the simple quadratic structure of the forward operator convergence 
rates results are difficult to obtain. Indeed, autoconvolution ‘scrambles’ the input 
function in such a way that qualihed nonlinearity conditions like the tangential cone 
condition and classical source conditions are generally not satished. For details of 
the corresponding results we refer to [m UHl IMl El] and the recent paper [7]. 
Moreover, we refer to mm for alternative approaches to deautoconvolution. 

A completely new approach to the autoconvolution problem from a mathemati¬ 
cal point of view was started in 2011 at the TU Chemnitz in collaboration with the 
research group ‘Solid State Light Sources’ of the Max Born Institute for Non-linear 
Optics and Short Pulse Spectroscopy, Berlin. The new onset was motivated by 
the development of the SD-SPIDER (self-diffraction spectral phase interferometry 
for direct electric-held reconstruction) technique in ultrashort (femtosecond) laser 
pulse characterization at the Max Born institute. A mathematical model for this 
method can be formulated by generalizing (|^ as a kernel-based autoconvolution 
equation 


( 2 ) 


k{s,r) f{s - t) f{T)dT = g{s), 0 < s < 2, 


with complex-valued functions / over M with supp(/) C [0,1] and g over M 
with supp(5f) C [0,2], and a complex-valued kernel k{s,T) with (s, r) G and 
supp(fc) C [0, 2] X [0,1]. In the sequel we use the polar coordinate representations 

(3) /(r) = a(r) and g{s) = b{s) 

for the searched-for function / and the right-hand side function g, respectively. 
For a detailed explanation of the physical background, i.e., of the SD-SPIDER 
approach in light of ([^ and concerning the availability of optical measurements 
for the amplitude functions a,b and the phase functions (p,^’ in (j^, we refer to 
the subsequent Section 

To be precise, we consider ([^ as a nonlinear operator equation in a Hilbert 
space setting. Our focus is on the spaces 

(4) X:= 4(0,1) and Y := Ll{0,2) 

of complex-valued square integrable functions with associated norms || • ||x, || • ||y 
and inner products (•, ■)x, (•, ■)¥, respectively. Then we write (|^ in the concise 
form 

(5) 


F{f) = g, fex, geY, 


VARIATIONAL DEAUTOCONVOLUTION AND PHASE RETRIEVAL 


3 


where /'^ G X denotes exact solutions to (|^ for given exact right-hand side g = 
G X, and where the forward operator F : X ^ Y, taking into account the 
specihc support intervals of pre-image and image functions, is dehned as 

[F(/)|(s) = |B(/)|(i,)e‘l*M‘> 

min(s,l) 

( 6 ) r 

:= J fc(s,r)/(s-r)/(r)dr, 0 < s < 2, 

max(5—1,0) 

for preimage function /(r) = a(r) 0 < r < 1. From (|^ it is seen that the 

support of the kernel function k is contained in the parallelogram 

(7) ^ := {(s, r) G [0,2] X [0,1] : 0 < r < 1, r < s < r + 1}. 

According to physical models (cf. |H [HI [15]) the kernel functions occurring in 
laser pulse characterization are continuous complex-valued functions on We 
will thus restrict our considerations to kernel functions, such that 


(8) supp(A;) C fc G and k{s, r) = k{s, s — t) for (s, r) G 

The last assumption in ([^ indeed holds without loss of generality as we have, for 
arbitrary kernel functions k{s,T) and all s G [0,2], the identity 


min(s,l) 

J Hsx)f{s 

max(s—1,0) 


min(s,l) 

r)/(r)dr= J 

max(s—1,0) 


k{s, t) -|- fc(s, s 
2 



T)f{T)dT. 


We mention that the structure ([^ is in particular satished if the kernel is generated 
by a complex-valued function k with supp(«:) C [0,1] and k G C[0, 1] such that 

(9) k{s,T) = k{t) — t) if (r, s) G and /c(s,r) = 0 otherwise. 

For some statements in Section we will be forced to focus on the special case 

(§, 

In this paper, we distinguish two different inverse problems in the context of the 
forward operator F from (|^ under the assumption (|^ for the kernel k, which are 
briefly presented in the following. 


1.1. The deautoconvolution problem. We call the inverse problem of identi¬ 
fying the complex-valued function / solving the operator equation (|^ deautocon¬ 
volution problem. In this context, it is assumed that only noisy data of the 
complex-valued function g"^ with 

( 10 ) \\g^-g^\\Y<S 

are available for a reasonably small noise level 5 > 0. 

This deautoconvolution problem was tackled in [15] in a direct manner, where 
theoretic consideration were only made for the trivial kernel k = 1. For the sta¬ 
ble approximate solution including a nontrivial kernel, an iterative regularization 
approach based on a variant of the Levenberg-Marquardt method was used there, 
which showed quite good results in numerical case studies with synthetic data. In 
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this paper, we extend assertions on properties of the forward operator F and on 
the solution of the operator equation ([^ to kernels k from (|^ or (|^. 

The real data situation in laser optics, see for details Section below, has an 
advantage and a disadvantage. Fortunately, noisy measurements for the amplitude 
function a in the solution / (see (|^) can be provided. These data were exploited 
in na to control a regularization parameter. However, measurements for the 
amplitude function b in g, which were required in combination with data for the 
phase ip, are unfortunately not sufficiently reliable in practice. As a result, the 
method in [15] failed for real data from optical experiments. In order to resolve 
this shortcoming, we clarify the objective and the data situation in the following 
manner. 

1.2. The phase retrieval problem. We call the inverse problem of identifying 
the phase function (p in the solution / (see ([^) of equation (|^ from noisy data of 
the phase function pj in the right-hand side g phase retrieval problem when noisy 
data of the amplitude function a in /, but no data of the amplitude function b in 
g are available. 

A hrst trial of a more adapted regularization approach with focus on the phase 
retrieval problem was included in the recent paper |3| taking into account the 
real data situation. However, the well-posedness of the non-standard variational 
regularization methods was not considered there. We are going to close this gap 
in the present paper, in addition to an essential improvement of the numerical 
implementation based on NURBS (Non-Uniform Rational B-splines; see, e.g., |24j ) 
in combination with a TIGRA-tjpe algorithm (the name is derived from Tlkhonov- 
GRAdient method; cf. [25l [26l l3^ ). 

The paper is organized as follows. After an explanation of the background 
from laser optics in Section we summarize properties of the forward operator 
F as well as of the deautoconvolution problem solving equation (|^ in Section 3 
Well-posedness results for the phase retrieval problem are presented in Section 4 
Turning to numerical considerations, we then briefly introduce planar NURBS anc 
describe how they may be utilized in the context of the complex-valued deauto¬ 
convolution problems in Section]^ Finally, in Section]^ we show numerical results 
obtained with a TIGRA-type method for both synthetic and real data. 

2. Physical Background 

The motivation for the complex-valued and kernel-based deconvolution prob¬ 
lem of solving (|^ in the introduction is the SD-SPIDER (Self-Diffraction Spectral 
Phase Interferometry for Direct Electric-held Reconstruction) method in laser op¬ 
tics. The aim of this method is the reconstruction of the electric held E{t) of 
ultrashort (femtosecond) laser pulses, which is a real-valued oscillatory function 
of time t. This function is usually decomposed into an amplitude part and an 
oscillating part by 

(11) 


E{t) = cos{a;ot + git)}, 
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where I{t) is the intensity (up to re-normalization), uo the carrier frequency, and 
ri{t) is called temporal phase. Since measurements are available for the spectral 
domain only, we have to consider the Fourier transform of E{t) 


OO 



— OO 


as a complex-valued function of the frequency oo. The Fourier transform can be 
written in polar coordinates as 

(13) ^(o;) = 

where S{u) is the spectral power density (up to re-normalization), and 0(a;) is 
called spectral phase. Fortunately, one is able to measure the spectrum directly, 
but contaminated with noise. Hence, the approximate identihcation of the desired 
physical quantity E{t) by inverse Fourier transform requires hrst the determination 
of 0. Careful denoising, e.g., by Fourier hltering or adjacent averaging of the 
measurements for S may prove helpful to reduce a possible oscillatory behavior of 
the algorithms. Since the spectral phase 0 cannot be measured directly, nonlinear 
optical processes have to be employed to infer the spectral phase from an indirect 
measurement. To this end, one suitable process is self-diffraction (SD), which is a 
spectrally degenerate variant of the more general four-wave mixing process. The 
electric held of the generated pulses in the spectral domain F^sd(iv) is related to 
E{uj) through 

(14) ^sd(w) = j K{oj,n)E{u + uc^-n)E{n)dn. 

0 

The SD process involves interaction with a continuous wave at a known frequency 
Ucvi, which has been incorporated into the kernel function K{uj,Q). Moreover, 
the kernel function K{u,Q) as an apparatus function is in principle known from 
physical modeling and can be assumed to be smooth in both arguments u and 
D. For a detailed description of the kernel see HI. We also write Esb in polar 
coordinates as 

(15) Esniu) = 

Perfect knowledge of all efficiency calibration factors in K provided, one can di¬ 
rectly use the spectral power density S'sd for reconstruction of E{t). Even a rela¬ 
tively small miscalibration, however, may have a dramatic inhuence on S'sd, which 
is why algorithms relying on the spectral power density may have trouble to pro¬ 
duce reasonable results [H]. The phase 0 sd, on the other hand, remains widely 
unaffected from such amplitude calibration problems. The measurement setup 
for the SD-interferogram is shown in Figure For physical details we refer to 

[laiisiEsi. 
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Figure 1. Measurement setup in self-diffraction spectral interferometry. 


The measurements for the spectrum of E indicate that the corresponding func¬ 
tion can be neglected outside a compact interval [wiow, ^^up], i.e. we assume that 

(16) supp E C [a;iow,Wup]. 

Consequently, Est> can also be neglected outside the compact interval 

(17) SUpp F/sD C [2cUiow ^cw) 2CiIup CJ^w]- 

Thus, (14) can be written as 

minfajup ,(.j-|-tJcw —t^iow } 

(18) ^ sd ( w ) = J K{u,n)E{u + uc^-n)E{n)dn. 

max-[tJjow —^up } 

The substitutions 

. ^sd( ^cw “t“ 2cUiow “b '^(^up ^low)): 

^('^7 * -^( ^cw “b 2cUiow “b '^(cJqp ^low)? ^low “b '7~(^up ^low)) 


reformulate (18) as 

min{l,s} 

(19) j k{s,T) f{s -r) f{T)ds = g{s), 0 < s < 2, 

max{0,s—1} 

which is precisely the mathematical model presented in the introduction. The 
physical background justifies the continuity of the kernel F in ^ as required in 

3. Some properties of complex autoconvolution operator and 

DEAUTOCONVOLUTION PROBLEM 


In this section, we summarize properties of the complex autoconvolution oper¬ 
ator and discuss the classical Tikhonov regularization for the deautoconvolution 
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problem aimed at solving the operator equation ([^ based on noisy data under 
the noise model (10). 


The Frechet differentiability of the autoconvolution operator mapping in the 
real Hilbert space L^(0,1) and the structure of its Frechet derivative were out¬ 
lined in [in]. Under the kernel assumption ([^ imposed on the kernel k, such 
result can also be formulated for the complex case § with the Frechet derivative 
F'(/) : X = 1) —y = 2) given for all / G X by the formula 

min(s,l) 

(20) [F\f)h]{s) = 2 j k{s,T) f{s — t) h{T) dr, 0 < s < 2, h G X, 

max(s—1,0) 

from which we easily derive 


(21) [F{f + h)- F{f) - F\f) h]{s) = [Fih)]{s), 0 < s < 2, h G X. 

Hence, we obtain with fp from ([^ and using 

(22) k ■.= max |fc(s,r)|, 

(s,T)e‘p 

the norm equation 

\\F{f + h)-F{f)-F'{f)h\\y = \\F{h)\\Y for all /, h G X 
and the nonlinearity condition 

(23) ||F(/ + h)- Fif) - F\f) h\\y < k \\h\\j, for all /, h G X. 

Such condition, which was already used in jO] Section 10.2] to obtain convergence 
rates for Tikhonov regularization of nonlinear operator equations, is by now the 
only available nonlinearity condition for the operator F from (|^. 

Note that the adjoint operator F'{f)* : Y —)■ X can also be written down 
explicitly, namely as 

T+1 

(24) [F'(/)V](r) = 2 j k{s,T) f{s - t) r{s) ds, 0 < r < 1, feX,reY. 


As usual we apply in the sequel the symbol B{f,r) for the closed ball in X with 
center / and radius r > 0 and the symbols ^ and —)■ for weak and norm conver¬ 
gence, respectively, in the occurring Hilbert spaces. The proofs of the following 
two results have been postponed to the Appendix. 

Proposition 1. Suppose that the kernel k satisfies Then the autoconvolu¬ 
tion operator F : X ^ Y from ^ is weakly continuous, i.e., for every sequence 
{fn}nm e X with fn^ fo in X as n ^ oo we have F{fn) X(/o) in Y. 

The following proposition outlines a very specific property of the autoconvo¬ 
lution operator, namely that the nonlinear F : Lp(0,1) —?• L^(0,2) from (|^ is 
non-compact, whereas the linear Frechet derivative operator F'{f) is compact in 
all points / G Lp(0, 1). This is a remarkable property of the autoconvolution op¬ 
erator, which rarely occurs for nonlinear operators. Conversely, it is well-known 
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that the Frechet derivative of a compact nonlinear operator is always compact. 
Again the proof can be fonnd in the Appendix. 


Proposition 2. Suppose that the kernel k ^ 0 satisfies (8). Then the kernel-based 
nonlinear autoconvolution operator F : X ^ Y from (6) is not compact. More 
precisely, there exists a sequence C 1) such that ^ 0 hut hn fit ^ 

in X as n ^ oo, which implies that we have, for all f & X and all r > 0, a 
sequence {/„ := f + rhn}nen C B{f,r) with fn^fin X, F{fn) F{f) in Y but 
^ifn) fit F{f) in Y as n ^ oo. On the other hand, the Frechet derivative F\f) 
given by (20) is a compact linear operator for all / G X. 


Remark 3. The proof of Proposition is based on finding seqnences ^ / in 
X snch that F{fn) fit F{f) in Y asn—)-oo. A conversely related concept is the 
local ill-posedness of the operator eqnation (|^ at / G X, which reqnires sequences 
fn in any neighborhood B{f,r) of / such that /„ fit f, but F{fn) — )■ F{f). For 
the specihc case fc = 1 on an example of the form 

fn = f + rhn, with h„(r) = 


was provided in [71 Example 3.2], showing that the operator equation ([^ is locally 
ill-posed everywhere in this situation. It is worth noting that, while the freedom 
in choosing r > 0 is exploited here for local arguments with arbitrarily small r, 
it may also be used for construction of elements /„ at arbitrarily large distances, 
11/ ~ fn\\x = 'f') whose images F{f) and Fi^ffi) are virtually indistinguishable, say 
||F(/) — F{fn)\\Y < This means that highly oscillating perturbations imposed 
on / cannot be recovered by a simple least-squares approach from noisy data 
oi g = F{f) with noise model (10), even if the noise level 5 > 0 is arbitrarily 
small. Since the same phenomenon has to be expected for more general kernels k, 
a regularization method seems to be always required in order to avoid oscillating 
numerical approximations. 


These observations are particularly interesting in combination with the unique¬ 
ness assertion in Proposition below which is based on the following well-known 
Titchmarsh convolution theorem; cf. |3T]. 

Lemma 4. Let ^ 1,^2 ^ with supp(^/) C [0, cx)), / = 1,2, and let for some 

constant a > 0 

S 

J ^i(s — r) ^ 2 ('r) dr = 0 for almost all s G [0, a]. 

0 

Then there are nonnegative constants Oi and 02 such that ai 02 > a and 

.^i(r) = 0 for almost all r G [0, oi] and ^ 2 (t) = 0 for almost all r G [ 0 , 02 ]. 

Proposition 5. Let the kernel k from he generated by a function k G ^[O, 1] 
with supp(fi:) = [0,1] such that 0 is satisfied. If for given g ^Y = L^(0,2) the 
function /I G X = Lp(0,l) solves then p and —p are the only solutions of 
this equation for the right-hand side g. 
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Proof. Based on Lemma the assertion of this proposition can be shown in 
analogy to the proof of Theorem 4.2 in [15], taking into account that under the 
stated assumptions on n the equality k{t) /^(t) = k{t) /(r) with / G X and for 
almost all r G [0, 1 ] implies that /I and / are the same elements of L^(0, 1 ) when 
taking into account that 

supp(k) = [ 0 , 1 ] holds if and only if k 7 ^ 0 almost everywhere on [ 0 , 1 ]. 

□ 


Note that Lemma under the assumption ([^ also provides us with a neces¬ 
sary and sufficient condition for the injectivity of the Frechet derivative; see (20). 
Namely, the bounded linear operator F'{f) : X —)■ F is non-injective if either / is 
the zero function almost everywhere on [0,1] or if k = 0. Otherwise the Frechet 
derivative is injective. 


Variational (cf. [28]) and iterative (cf. [22]) regularization methods form two 
classes of standard methods for stabilizing ill-posed nonlinear operator equations 
(|^ in Hilbert spaces. The most prominent representative for the hrst class is 
the Tikhonov regularization, where in the simplest form (cf., e.g., [Q] Chapt.lO]) 
regularized solutions are minimizers of the extremal problem 

(25) ||F(/) - /||y -1-a 11/-/III-)■ min, subject to / G X, 


with some initial guess / G X. For obtaining convergence rates of the regularized 
solutions, an appropriate interplay of solution smoothness and structural condi¬ 
tions expressing the nonlinearity of F in a neighborhood of the solution is required 
(see, e.g., m for an overview). For the Tikhonov regularization of the form ( [2^ 
and F from the condition (23) acts as nonlinearity condition sufficiently well 


and allows proving the convergence rate 


\\fiiS)-f\\x = 0(y6'^ as 5^0 and «(5) ~ 5, 

since the operator F is weakly continuous (cf. Proposition [^, which implies that 
F is weakly closed. The latter result additionally requires the existence of a source 
element w G Y satisfying the smallness condition ||tc||y < 1 such that the source 
condition 


T+1 


(26) 


/*(’■)-/(’■)= Hs.t) f(s - T)w{s)ds, 0<T<1, 


is fulhlled (cf. [HI Theorem 10.4]). However, in [71 Proposition 2.6] it was shown 
that such source condition (26) is hardly possible to achieve even in the simplest 
cases of a kernel k = 1. 


For wide classes of iterative regularization methods, however, the tangential 
cone condition 

(27) \\F{f + h) - FU) - F'(f) AIK- < c||F(/ + h) - F(/)||k tor all h e B(f,r) 
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is required, with a constant 0 < c < 1 , some radius r > 0 , and at least for / G X 
in a neighborhood of the solution p to (|^ (see, e.g., P Chapt.ll]). Slightly 
modihed variants of the nonlinearity condition (27) for constants 0 < c < cxd and 
terms 9{\\F{f + h) — F(/)||y) instead of \\F{f + h) — F(/)||y, with concave strictly 
increasing functions 6 : ( 0 , cx)) — )■ ( 0 , cxd) and lim 6{t) = 0 , are also relevant for 

obtaining convergence rates in Tikhonov regularization, in particular if the solution 
smoothness is low and approximate source conditions apply (see m Section 4.2]). 

Proposition 6 . For the autoconvolution operator F : X ^ Y from the 
tangential cone condition (21) cannot hold for any f E X with a constant 0 < c < 1 
and a radius r > 0. 


Proof. The assertion is, in principle, a consequence of the noncompactness of F 
in contrast to the compactness of F'{f) for all f E X. Using the triangle inequality 

||F(/ + h)- F{f)\\y < \\Fif + h)- F{f) - Fff) h\\Y + \\F\f) h\\y, 


we would have from (27), for hxed f E X and r > 0, the inequality 

(1 - c)||F(/ + h)- F{f)\\y < \\F\f) h\\y for all h E S(/,r). 


By substituting h := hn in that inequality with C I3(f,r) from Proposi¬ 

tion [^satisfying hn ^ 0, hn -p 0 in X and F(f + hn) -p F{f) inY as n —)■ cx), we 
arrive at a contradiction, because F'{f) is compact and hence satishes the limit 
condition lim \\F'{f) hn\\y = 0. This contradiction proves the proposition. □ 


Remark 7. As originally discussed in [19], the four norm terms associated with 
the Taylor remainder equation, 

F{f + h)- F{f) = F'if) h + Rem(/, h). 


namely ||F(/ + h) - F(/)||y, ||Rem(/,h)i|y = ||F(/ + h) - F{f) - F'{f)h\\y, 
||F'(/) h||y, and ||h||x show distinguished cross connections depending on whether 
the corresponding nonlinear operator equation (j^ is well-posed or ill-posed. The 
convergence of well-posed problems usually results from the fact that the remain¬ 
der ||Rem(/, h)||y converges to zero faster than the term \\F{f + h) — F(/)||y as 
ll^ilx —t 0. In the ill-posed situation, however, where F is a ‘smoothing’ operator, 
the term ||F(/ + h) — F(/)||y may be signihcantly smaller than ||Rem(/, h)||y 
even for arbitrarily small ||h||x- In the latter case, there exist operators F for 


|x- 

which the tangential cone condition (27) may fail to hold even for large constants 
c > 1 as well as for modihcations with terms 6{\\F{f + h) — F(/)||y) rather than 
||F(/ + h)-F(/)||y. Then only Lipschitz continuity conditions for F'{f) such 
as (23) can be seen as nonlinearity conditions. In this particular situation, the 
linear operator F'{f) does not cover sufficient information about the nonlinear 
operator in a neighborhood of B{f, r) for ensuring convergence rates of regularized 
solutions. This information dehcit seems to be the case for the autoconvolution 
operator F from (j^ and, unfortunately, prevents the establishment of convergence 
rates in Tikhonov regularization when the standard source condition (26) fails. 
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4. A REGULARIZATION APPROACH FOR THE PHASE RETRIEVAL PROBLEM 


In this section, we collect some basic well-posedness results concerning the regu¬ 
larization approach by means of minimizing Tikhonov-type variational functionals 


for the different problems summarized in Subsections 1.1 and 1.2 


The theory of Tikhonov-regularization is well understood for ill-posed problems 
of the general form 


(28) 


J^(x) = y, 


where a; is a searched-for quantity of interest and y is approximately known from 
measurements in the form of data y^. We refer to the monographs [siEsiEni as 
well as to the seminal works m [T8] for detailed regularization results. For our 
purposes, let us recall the main sufficient condition on T : dom(J^) d X [2HI 
Section 3.2] for the well-posedness of minimizing a Tikhonov-type functional 

'Tiix) = ^|l-F(a:) - + aR{x) 

in Hilbert spaces X,y: namely, that the forward operator X be weakly sequentially 
closed in the sense that 


Xn ^ X & X and X{xn) ^ y E y x E dom(J^) and X^x) = y 

hold for all sequences {x„} C dom(J^). 


Remark 8. Note that we do not specify further the choice of the penalty term 
7l{x). For the following results to hold true, 7l{x) is required to be proper, weakly 
sequentially lower semicontinuous and to have weakly sequentially precompact 
sublevelsets. We will summarize these properties by saying that 7l{x) is assumed 
to be stabilizing. The interested reader will easily verify that for stabilizing penalty 
terms TZ{x) in combination with weakly sequentially closed forward operator X 
with dom(J^) = A, Assumption 3.13 in [28] is satisfied (with the exception of 
convexity of 'R. which is, however, not required here). Consequently, the well- 
posedness results in [281 Section 3.2] hold true. In particular, we have: 

• Existence of minimizers (cf. Theorem 3.22); 

• Weak, subsequential stability of the minimizers (cf. Theorem 3.23); 

• Weak, subsequential convergence of the minimizers to an 7^-minimizing 

solution as ^ 0 under suitable parameter choice rules (cf. Theorem 3.26). 


The deautoconvolution problem of Subsection O is naturally of the form ( [2^ 
with T = F from (§. and due to the properties in Section we immediately 
obtain the following regularization result. 


Proposition 9. Suppose that the kernel k satisfies and that the penalty term 
Tl{f) is stabilizing in X = ^^(0,1). Then for any g°E V = T^(0, 2), minimizing 
the Tikhonov functional 

for f E X is well-posed in the sense of Remark\^ 
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Proof. The weak sequential closedness of F : X ^ Y follows readily from the 
weak sequential continuity in Proposition [Tl The result then follows arguing as in 
Remark I □ 


In the particular situation of phase retrieval problems as in Subsection 1.2, the 
data consists of two separate parts. On the one hand contains information 
about the absolute value of the unknown function /, and, on the other hand, 
'0'^ arg(i/l) is related to the image space Y. For formulating this problem in the 
context of Tikhonov regularization, we first dehne the Sign operator on L^(0,2) 
pointwise almost everywhere (a.e.) by 


9{s) 

[Sign(c/)](s) = <j 


if g{s) ^ 0 
else. 


Then the forward operator in the sense of (28) also consists of two parts and maps 
as 


J-PR : X = Lj(0,1) ^ = Lj(0,2) X L"(0,1) 

TMf) = (Sign(F(/)), I/I), 

where the function |/| G L^( 0 , 1 ) is again defined pointwise a.e.by |/|(r) = |/(r)| 
for f E X. Note that the product space y = YxL‘^{0, 1) is a Hilbert space endowed 
with any of the equivalent norms ||(2/, := Ibll Y + /3||®IIa corresponding to 

discrepancy terms 

||Fpr(/) - ||Sigti(F(/)) -e‘*‘||?. + /3|| |/| - V|li. 

Even though the latter functional by itself resembles a Tikhonov functional, we 
emphasize that here the factor f3 > 0 acts as a balancing weight between two 
discrepancy terms and does not assume the role of a regularization parameter. This 
difference is reflected in our notation, where in the following we clearly distinguish 
between regularization parameters (denoted by a) and the discrepancy weight 
/9. Nevertheless, it proved reasonable to consider well-established regularization 
parameter choice rules also for choosing f3, and we present a suitable example in 
Section [6l 

The operator : X —)■ 3^ is, however, not weakly sequentially closed as 
neither Sign( 5 f) nor |/| have this property, which is easily seen by considering the 
sequences Qn = \ and fnif) = e®"'*, respectively. Thus, following the approach 
suggested in [T3], we approximate the Sign operator for £ —)■ -|-0 by 

While the latter operators are Lipschitz-continuous with constant they still do 
not satisfy the sufficient condition for well-posedness of Tikhonov regularization. 

Lemma 10. For any £ > 0, the operator Signj,( 5 f) : L^(0,2) —)■ L^(0,2) is not 
weakly sequentially closed. 
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Proof. To prove the assertion we construct a counterexample. Let K > e he 
fixed and define the sequence {gn} recursively by 


9o{t) 


2K if 0 < t < 2/3 
-K if 2/3 < t < 2 


and gn+iit) 


gn{2,t) if 0 < f < 1 

gn{2t — 1) if 1 < f < 2. 


In fact, for these choices gn ^ 0 holds, yet Signg( 5 f„) ^ | 7^ Signg(O), which 

shows that Sign^ is not weakly sequentially closed. □ 


In a sense, the weak topology in L^(0,2) is insufficient to derive continuity 
results even for Sign^( 5 '). A common remedy for obtaining regularization properties 
for operators that are continuous is to restrict ourselves to some subspace with 
compact embedding into Lp(0,2). For the phase retrieval problem, we may thus 
formulate the following well-posedness result. 


Proposition 11. Suppose that the kernel k satisfies Q) and that the penalty term 
7?.(/) is stabilizing in Xq = Then for any e > 0, E L^(0,1) and 

G L^(0,2), minimizing the Tikhonov functional 


(29) TfijU) = 2 Sign,(F(/)) 

for f E Xq is well-posed in the sense of Remarkl^ 


sr 


2 d II 


II „ + aTZ(f) 


Proof Due to the compact embedding of Xq X, both operators 

Sign,(F(/)) : Xo ^ y and |/| : Xq ^ X 

are strongly sequentially continuous. Hence they are in particular weakly sequen¬ 
tially closed. The result then follows from the identical argument as in Remark 

□ 


Remark 12. As our notation suggests, the measurement errors affect the phase 
function data additively, i.e., 

^jJ^ = ^jJ^ -\-1), 

where -d E L^(0,2) denotes the noise. Due to discontinuities in the principal part 
of the complex argument, however, we evaluate data discrepancy on the complex 
unit sphere S^, where we have a multiplicative noise model, 

giy* _ g*y+. g*i? 


Finding a discrepancy term that suitably addresses this particular situation is cer¬ 


tainly future work, but motivated by the regularization results in Proposition 11 


and the successful numerical experiments in Section we use the subspace topol¬ 
ogy from Y = Tc(0, 2) on : -0 G L^(0, 2)} C X instead, and hence introduce 
the pseudo-metric d{fii,'ijj 2 ) ■= for '0i,'02 G L^(0,2). Note that this 

distance penalizes phase differences modulo 2n and does not increase the noise 
level, as 


d{fi ,0T) = 


2 sin 


L2(0,2) 


< 110 "^ - 0 ^ 


Il2(0,2)' 
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As mentioned in Section the laser pulses generated in optical experiments 
exhibit a limited bandwidth; cf. (16). In the context of our current notation 
this requires that the searched-for function /I (which corresponds to the Fourier 
transform E of the laser pulse) can be neglected outside a certain parameter range 
[ti, C (0,1). Due to the structure of the autoconvolution operator, the same is 
to be expected for the image = F{p). Therefore, the phase data ~ arg(i/l) 
carries little or no useful information in regions where |i/^| is close to zero. As 
reliable measurements for |i/l| are, however, not at hand, this fact is not accounted 


for in the Tikhonov-functional (29). To overcome this issue, a different data fidelity 
term has been proposed in |1], where the phase data discrepancy is weighted by 
|F(/)| or, more precisely, by the normalized function |F(/)|/||F(/)||y. While this 
normalization is necessary to avoid an artificial bias towards reconstructions for 
which |F"(/)| is small, it also introduces a singularity for ||F"(/)||y = 0. Therefore, 
we again introduce an approximation level £ > 0 and consider the variational 
functional 
(30) 


vUf) = 






WFimi 


+ 2 II I/I ~ 


+ 00 , 


l^ + aW), ii\\Fif)\\y>e, 
else. 


with Y = L‘^{0,2). The proof that this functional also admits a minimizer in 
H^{0, 1) is included in the Appendix. 

Proposition 13. Suppose that the kernel k satisfies ^ and that the penalty term 
is stabilizing in Xq = in the sense of Remark^ Then for any 

a®" G T^(0, 1), & -f"^(0, 2), and £ > 0 such that 


dom7^(/)n{/eXo : ||F(/)||y > 4 ^ 0, 


a minimizer of ^{f) defined by (30) exists in Xq. 


Remark 14. Another phenomenon of ill-posedness in solving the autoconvolu¬ 
tion equation (cf. [HI [15] and [71 Example 3.2]) is due to amplitudes a(r) in 
/(r) = air) occurs when a(r) blows up to infinity locally near some 

To G [0,1] in a way that hardly leaves effects on F{f). This phenomenon, how¬ 
ever, is suppressed in case of the phase retrieval problem because the amplitude 
function is known up to measurement errors. It was, indeed, proven in [71 
Proposition 3] that even locally well-posed situations arise under the assumption 
that the amplitude function a be fixed and essentially bounded. 


5. Non-uniform rational B-splines (NURBS) 

The numerical solution of the deautoconvolution problem |1.1 as well as of 
phase retrieval problems L2 requires discretization of the complex-valued auto¬ 
convolution equation with kernel function (|^. The most natural choice are dis¬ 
cretizations using piecewise constant functions, either in terms of step functions 
{X[i/n,(i+i)/n)}i=o,...,n-i (scc, c.g., [miHlElH] ) or by means of Haar wavelets (e.g.. 
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HI [32]). While step functions yield simple (and computationally efficient) formu¬ 
lae, essentially reducing the continuous autoconvolution to its discrete counter¬ 
part, Haar wavelets are particularly suitable for the reconstruction of functions 
in L^( 0 , 1 ) as they yield an orthonormal basis both of the inhnite dimensional 
Lebesgue space as well as of its truncated, hnite-dimensional approximations. Mo¬ 
tivated by the results in Section concerning the existence of minimizers of the 
Tikhonov functional in the smoother space Tfc(0,1) and by the underlying physi¬ 
cal problem of ultrashort laser pulse characterization (see Section]^, we focus on a 
different representation in the sequel, which is taylored towards the reconstruction 
of smooth functions. Namely, we will model the curve ( Re /, Im /) in the complex 
plane (and hence / itself) as a rational B-spline curve. 

The shape of a non-uniform rational B-spline (short: NURBS) is determined 
by a set of n control points {Pi,...,P„} C corresponding positive weights 
w = {tCi, ... ,Wn} and a non-decreasing knot vector rj = {rji, ..., rjn+p+i}, where p 
is the polynomial degree of the spline. From the knot vector i], the B-spline basis 
functions Nj^p are determined via the Cox-de Boor recursion 


/V iri - / ^ if hi < ^ < hi+i 
1 0 otherwise. 


NjAr) = 


T — 


Vj 




Vj+p+1 


'Nj+l,p-l{'T)- 


Vj+p Vj Vj+p+i Vj+i 

Notice that throughout this section we adopt the convention 2 := 0 as is customary 
in this context. The NURBS curve is then given by 

Ei=l PjWjNj,p{T) 


7 [P,u;](r) = 


Y.U^,N,p{r) ’ 


Lmin ^ T ^ Rnax) 


where Tmin := hp+i, Umax := Vn+ 1 - We refer the interested reader to [211 EZ] for a 
comprehensive introduction to NURBS. Defining rational basis functions as 

WjN,^p{t) 

the curve can be equivalently written as 

n 

v[P,w]{t) = ^P,Pj-p(r). 
i=i 


It can be easily seen from these definitions that NURBS curves are invariant (up 
to re-parametrization) under rescalings and shifts of the knot vector. Without 
loss of generality, we may thus assume that the knot vector satishes rjp+i = 0 and 
rjn+i = 1, so that the NURBS curve is parametrized by r G [0,1]. A common 
choice are open knot vectors of the form 


V = (0,...,0,7p+2,...,hn,l,---,l), 

which have the additional property that the resulting NURBS curve 7 (r) begins 
and ends in the first and last control point, respectively, i.e., it satisfies 7 ( 0 ) = Pi 
and 7 ( 1 ) = Pn- 
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Autoconvolution and NURBS. For planar curves the NURBS control points 
Pj are typically assumed to be in Nevertheless, the formulae and results in 
Section remain valid, if we model them as Pj G C which allows to conveniently 
view the parametrized spline curves as complex-valued functions. In addition, 
we shall refer to the real and imaginary parts of the control points by Uj and 
Vj, respectively, so that Pj = Uj -\- ivj. The main motivation for working with 
real parameters is the fact that real-valued functionals, such as the variational 
Tikhonov-type objective functionals dehned in Section are not holomorphic. 
When using gradient based optimization methods, it is therefore necessary to 
consider partial derivatives with respect to u and v. 

Now let the spline degree p, the number n > p -|- 1 of control points, as well 
as the knot vector rj be hxed. Then, we denote the hnite dimensional space of 
NURBS design parameters x = {u,v,w) by 

Xn := M” X X C 

and by 7 : —)■ X = L^{0, 1 ) the synthesis operator 


7W(^) = + iVj)RjA'r), 


with RjAA 




In terms of the complex-valued NURBS curves, the autoconvolution operator maps 
as 


(31) Fn — F o j : Xn — )■ U — Tc(0, 2), 

where F : X —)■ U is given by (§ with / = 7 [a;]. In this way, we aim to reconstruct 
approximations / in the discretized space 

X7:={7[x] : l6X„}cffi(0,l). 


Regularization for NURBS design parameters. As illustrated in Section 
the regularity of the forward operator for the phase retrieval problem 1.2 in its 
continuous form guarantees existence of minimizers of Tikhonov-type functionals 
in spaces that embed compactly into !)• choosing a discretization space 

that consists of sufficiently regular spline curves, we readily ensure that the re¬ 
constructed solutions belong to Rc(0,1). As penalty term for the NURBS design 
parameters x = {u,v,w), dehned in the previous subsection, we propose 


(32) 


7l{x) ^ j3pTlp(u,v) 0wTZw^(w), 


where we penalize the distance between control points Pj = Uj + ivj, and hence to 
some extent the length of the curve (corresponding to the 1 ) seminorm), by 




n—1 


^ ffi (u+i u) 


-t=i 
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To account for the constraint that the NURBS weights are required to be positive 
we penalize them by 


T^woiw) — ^ /, 


Wo 




i=i 


where 


{ +00 if tc < 0 

w~^ if 0 < w < wo/2 

{wo/2)~^{w — wo^ ifw>wo/2. 

Clearly, this functional, which is shown in Figure]^ acts as a (quadratic) barrier 
for w —+0 as well as for w —)■ +oo. In addition, it gives preference to values 
near a reference weight tco > 0. It is worth noting, however, that wo does not 
represent a priori knowledge on the average size of the weights, because NURBS 
weights are only determined up to a constant factor. Instead, the choice of wq 
in combination with the size of the parameter typically impacts the relative 
difference max(t(;)/ min(t(;) in the approximate solutions. 



Figure 2. Regularizing functional /«,o(u') for wq = 10 (left) and its 
log-log plot (right). 


6. Numerical results 


We have tested our method both on synthetic data as well as on real data ob¬ 
tained from an SD-SPIDER apparatus at the Max-Born-Institute for Nonlinear 
Optics and Short Pulse Spectroscopy in Berlin, Germany. The regularized solu¬ 
tions are NURBS design matrices = (Wq^, minimizing a variational 

Tikhonov-type functional 


' a,P 


(x) = 


F„(x) - |Fn(x)|e 




Y 




+ a|| ItMI -a'"||^ + a7^(x), 


with penalty term given by (32). The corresponding functions in X = Lf- (0, 1) 
are obtained as = 7 [x^’^]. To be precise, we used the functional (30) with 
approximation level e = 10“^°, but as situations with ||F„(x)||y < e never occured 
during the computations, we shall omit this additional parameter for the sake of 
brevity. 
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For finding an approximation of we use a Quasi-Newton method, 
(33) Xj+i = Xj + {xj), 


with step-sizes pk satisfying the Wolfe conditions and BFGS-updates of the ap¬ 
proximation Hk of the Hessian. It is well-known, that gradient based optimization 
methods for nonlinear problems generally suffer from local minima To obtain 
a better approximation of the global optimizer, we employ a strategy known as 
TIGRA (the name being derived from TIkhonov-GRAdient method), which was 
introduced by Ramlau ESI. This method was proven to converge globally for 
suitable problems and with suitably chosen parameters in [25l [26l |32] . 

Aiming at a well-balanced data £t both for ^ arg(i/l) and for ^ |/1| = 
| 7 [xl]|, we used a similar approach to TIGRA, but with respect to the weight 
balancing the discrepancy terms (cmp. Section]^. Starting with large /^o and 
g < 1 (here: /3o = 100, q = 0.25), the resulting method is described as follows. 


Find G arg min (x) using (33) starting from Xfc,o = 


< 7,5 




where (3k = q(3k-i- Initially, for large (3q only is emphasized while reconstruc¬ 
tions will not typically provide a good match for . Then the weight is gradually 
shifted giving more and more importance to while the initial good match for 
declines only mildly. A globalization approach (such as the proposed reweighing of 
the discrepancy terms, for example) is certainly required to reach an approximate 
solution of the original problem, but comes at the cost of a higher computational 
effort. In our case the additional effort is due to the repeated optimization with 
different values of the weight (3k- To keep the number of iterations to a minimum, 
however, we may solve the earlier optimization problems inexactly and increase 
the required precision while the iteration proceeds. In our experiments we have 
employed this technique using tol^ = max( 2 ^, 10 “®) as stopping tolerance for 
VT^’^^(xfcj), limiting the number of iterations during each optimization proce¬ 
dure to maxit = 10000. 

The opposing trends of the two data hdelity terms (cmp. Figure during the 
iteration can be exploited to obtain a stopping rule for the procedure. Gonsidering 
that our main objective is to optimize the overall data fit, we define a weighted 
relative least-squares functional 


(34) 


e(x)^ := 2 d{xY + r(x)^. 


in terms of the auxilliary quantities 


d(x)^ := 


Fn{x) - |F„(x)|e 




\Fn{.x) 


and 


r(x)^ := 


|7[x]| 


X — a 


|2 

lx 


lY 


12 

lx 


and stop the iteration when e(x^’^^) reaches its lowest value. Putting more em¬ 
phasis on the discrepancy term d{x) in Y proved benehcial in the numerical exper¬ 
iments. We remark that this procedure may be regarded as a bilevel optimization 
approach for choosing the discrepancy weight. Namely, for hxed a > 0, we would 
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Figure 3. Relative data error quantities d{x2^i^) (red) and 
(blue) from the following subsections in dependence on (3. The hnal 
value (3* is chosen as minimizer of e(x)^’^)^ (yellow). 



Figure 4. The real (left) and imaginary (right) part of the kernel function. 


choose 13* as an approximate solution of 

/ ( 7,(5 \2 

/3>0 

such that G argmin7)(^^(x). 

X&Xn 

For our purposes, however, the primary motivation is to improve the global con¬ 
vergence properties of the Quasi-Newton method rather than to solve the latter 
problem. The hnal values from the experiments in the following subsections have 
been collected in Table [U 

We hrst present some results obtained with synthetic data. In order to stay as 
close to the realistic data situations as possible, we use the same kernel function 
throughout the following subsections. The kernel shown in Figure]^ was obtained 
from physical modelling of the nonlinear optical processes that result in the mea¬ 
surements (cf. |23l [151 0]). The target phase function (pi = arg(/l) was chosen 
identical to m Section 3.2] (see also [I5]), and the noisy data were generated by 
adding 1% relative noise. 

The approximate solutions were represented as complex-valued NURBS curves 
/ = 7 [x] as introduced in Section]^ In our experiments we used quadratic splines, 
i.e., p = 2, with n = 150 control points and corresponding positive weights. The 
knot vector p, which determines where and how the control points impact the 
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parametrized curve, was chosen as open uniform. This is to say that 


Pi ■ ■ ■ Pp+l 0 ) Pn+l ■ ■ ■ Pn+p+1 I 5 


and 


Vj 


J -P 


n — p 


for p + 2 < j < n. 


In X = T^(0,1) we have discretized the resulting curves by choosing N = 1000 
equidistant sampling points and, correspondingly, in the image space 

Y = T^(0, 2) by Sm = for m = 1,..., 2N — 1. Computations were started 
with an initial guess consisting of constant weights Wj = Wq = 10 and control 
points Pj = Uj + ivj interpolated from a®" such that Uj = Vj ~ a/ 2 ■ In 


the penalty term (32) we also used wq = 10 and balanced both terms equally by 
l3p = l3w = 1. The regularization parameter was chosen as a = 10“®. 


Data 

CPU-time 

#Iter 

13* 




{a^,y^) 

1585 s 

19942 

3.81 ■ 10“^ 

1.02 ■ 10-2 

1.55-10-2 

4.47-10-^ 

(aUV'") 

1663 s 

19391 

9.54 ■ 10-® 

1.39-10-2 

1.68-10-2 

6.69 - 10-^ 

measured 

1848 s 

24403 

3.81 ■ 10-^ 

1.77-10-2 

1.97-10-2 

1.01 - 10-3 


Table 1. Comparison of results for various data situations from 
the following subsections with /3o = 100, a = 10“®. 


Phase retrieval with synthetic data for y. We first consider the ideal data 
situation, assuming that measurements are available for both the modulus and 
the argument oi ^ y as well as for the modulus a'^ ~ |// = |7[a:^]|. To 
simulate the real data situation we used the same measurements for as in the 
final subsection. These measurements are shown in Figure together with the 
reconstructed solution. 

In order to be able to work with the same parameter values for /q, tol^, and a 
as in the other test cases, we consider here the Tikhonov functional given by 

II ^ / X AI|2 


q-0-,d 


(x) = 


X 


y 


lY 


WyTy 


+ AII IbWI 


1^ + aTZ{x), 

where the NURBS penalty term is as defined by (32) and a = 10“®. Similarly, the 
relative data misfit term 


d{x) = 


\Fn{x) - y^ 


\Y 


\\y^ 


Y 


was used in the least-squares error functional e(x) defined by (34). 

The results are shown in FiguresJ^and]^ Note that the real and imaginary parts 
of the reconstruction //’^* = 7 [x^’^*] (right column of Figure 5) evidently provide 
a good approximation of /h Small oscillations in regions where the modulus |// is 
close to zero, however, result in quite large deviations of their arguments in these 
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Figure 5. Top-left: Target phase (red) and reconstruction 
arg(/^’^*) (blue). Bottom-left: Data a'^ |/'^| (red) and recon¬ 
struction (blue). Right column: Real (top) and imaginary 

(bottom) part of /'*' (red) and (blue) together with NURBS 

control points. 



Figure 6. Top: True arg(?/l) (red) and reconstructed argF„(x^’^*) 
(blue). The brightness is chosen in (logarithmic) dependence on 
|-^n(a;«(^.)l- Bottom: True \y^ (red) and |Rn(a;«(^.)l (blue) 


areas (left column of Figure]^. As one might expect reconstructions can therefore 
only be reliable in those regions where \ fap* \ 3> 0. To emphasize this observation 
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Figure 7. Top-left: Target phase (red) and reconstruction 
arg(/^’^*) (blue). Bottom-left: Data a®" |/1| (red) and recon¬ 

struction 1/^’^.I (blue). Right column: Real (top) and imaginary 
(bottom) part of together with NURBS control points. 


we have weighted the brightness of the plot of arg and arg(/l) depending 

on their absolute values. 


Phase retrieval with synthetic data for arg(y). As another academic ex¬ 
ample, we consider the phase retrieval problem corresponding to the real data 
situation in the following subsection. We assume that measurements are available 
for the argument oi ^ y and for the modulus ~ |/1| = | 7 [xl]|. Again, we 
used the available experimental data for a'^ which is shown in Figure 

The results are shown in Figures and When comparing to the ideal data 
situation where both |i/l| and |/1| are approximately known, a certain fall-off in 
quality is evident. However, especially in those regions where |/1| 3> 0, the recon¬ 
struction still provides a good approximation of the target. 


Phase retrieval with real data. In the real data from optical measurements, 
the frequency band containing the support of the a'^ ~ |/1| is located in between 
Wiow = 3.5 X 10 ^^ Hz and cUup = 4.1 x 10 ^® Hz, and the frequency of the continuous 
wave (see Section at cUcw ~ 3.86 x 10^^ Hz. 

As above, the Tikhonov functional is given by 

2 




(x) = 


Fn{x) - \Fn{x)\e 




/ 5 | 


\Fn{x) 




- + ||| I7NI 


a'^ll^ -|- aTZ{x), 


with penalty term 77(x) defined by (32) and a = 10 For the TIGRA-type 
approach with respect to the weight [3, we used /3o = 10^ and q = 0.25. Keeping 
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Figure 8. Top: Generated data (red) and reconstruction (blue) 
for argF„(a;^’^*). The brightness is chosen in (logarithmic) depen¬ 
dence on |F)i(x^’^.)|. Bottom: Target ||/^| and reconstructed 

l^n«’S.)|- 


in mind that phase data is reliable only in regions where the absolute value of the 
function is sufficiently large, we again observe a good data fit in Figures and 
To highlight this observation, we compare several reconstructions, which were 


10 


obtained using different parameters in Figure 11 


Conclusion 

In this paper, we have studied complex-valued autoconvolution problems with 
continuous kernel-functions in different data situations, arising, for example, in the 
characterization of ultrashort laser pulses by means of the SD-SPIDER method. 
We have derived fundamental analytical properties, in particular, weak-to-weak 
continuity in L^-spaces of the autoconvolution operator, which ensure well-posed- 
ness of regularization approaches by minimizing Tikhonov-type functionals either 
in Lp or in compactly embedded subspaces, depending on the data at hand. In¬ 
spired by the TIGRA method and using discretizations in terms of NURBS curves, 
we have also proposed a novel globalized numerical method for phase retrieval prob¬ 
lems corresponding to the real-world data situation, which is able to hnd solutions 
that appropriately reproduce the given data. 

Nevertheless, a number of open questions remain. Convergence rates results are, 
to be best of our knowledge, completely unavailable except for very special cases, 
and we have proven that classical nonlinearity conditions such as the tangential 
cone condition are not suitable to tackle autoconvolution problems. Concerning 
the noise model, further improvements might be possible with discrepancy terms 
that capture more adequately the multiplicative noise structure on the complex 
unit sphere for measurements of the complex phase function. Finally, we expect 















24 


S.W. ANZENGRUBER, S. BURGER, B. HOFMANN, AND G. STEINMEYER 



Figure 9. Top-left: Reconstructed phase arg(/^’^*). Bottom-left: 
Data a'^ ~ l/^l (red) and reconstruction (blue). Right colomn: 

Real (top) and imaginary (bottom) part of together with 

NURBS control points. 



Figure 10. SD-phase data (red) and reconstruction (blue) for 
aigFn{x^a*)- The brightness is chosen in (logarithmic) dependence 
on |F„(x^.)|. 


that an even better data £t could be achieved numerically by devising a method 
for blind deautoconvolution, i.e., by including the kernel function as a free or 
parameter-dependent variable in the optimization method. 
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Figure 11. Comparison of various phase reconstructions arg(/^’^*) 
for a = 10“® (blue), 5-10“® (yellow), and 10“® (red), showing a close 
match in the region where |/^’^*| is large. 
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Appendix 

Proof of Proposition We consider a sequence {fn}n&N ^ ^ with /„ ^ /o in 
A as n — )■ cx) and show that lim {F{fn) — F{fQ),ri)Y = 0 for all rj E Y, which 

n^oo 

proves the proposition. Indeed, by Fubini’s theorem we have 


{F{U)-F{fo),v)Y 


2 min(s,l) 

= J j k{s,T){fn{s -r)fn{r) - fo{s-T)fo{T))dTr]{s)ds 

0 max(s—1,0) 


min(s,l) 


k{s, T){fn{s - r) - /o(s - r))(/n(r) + fair)) dr 


max(s—1,0) 


7]{s) ds 


= J ifnir) + /o(r)) 
0 

1 

= j ifnir) + /o(r)) 
0 




Hs,T){fn{s - r) - /o(s 


T))ri{s) ds 


dr 


J k{^ + r, r)(/„(0 - fo{0)v{^ + 
.0 


dr. 


If we use the settings A„(^) := /„(^) - /o(0, 0r(O := + A'r)i7(^ + r) and 

Sn(T) := /o wliere 0^ G A for all r G [0,1] due to the continuity 
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of the kernel k and A„, G X for all n G N, we arrive at 

1 

{F{fn) - F{fo),ri)Y = j (/„(r) +/o(r))S„(r) dr = (/n + /o,S„)x. 

0 

By some calculations it can be shown that the family of functions n G N, 
is equicontinuous on the interval [0,1]. On that interval, the sequence 
converges pointwise to zero, because ^ 0 in A implies 

Sn('r) = (©r, A„)x —t 0 as n ^ oo for all r G [0,1]. 

However, an equicontinuous and pointwise convergent sequence of functions is even 
uniformly convergent, which yields lim^^oo ll^nllA = 0. Then the limit condition 

\{F{fn) - F{fo),r])Y\ = \{fn + fo,Fn)x\ < \\fn + foWxWFnWx ^ 0 aS 71 ^ OO, 

which is based on the fact that ||/n + /o||a is bounded, completes the proof. □ 


Proof of Proposition]^ In the special case fc = 1 on ip from ([^, the proof would 
be based on the fact that we have F{hn) 0 for := e*”* ^ 0. For a general 
kernel function A;, however, we have to take into account some more details. Since k 
is not identically zero, there is some (sq, Tq) G int(^) with k := A;(so, tq) 7^ 0. Due 
to the continuity of the kernel function /c on ^ there exists an open neighborhood 
Uq of (so,ro) with |/c(s,r) — A;| < -^ for all (s, r) G Wq © int(^). If we dehne 
D := {(2r, r) : 0 < r < 1}, then there exists (si, Ti) and e > 0 such that with 

Ui := {(s,r) G int(^) : Ei < r < Ti + £, Si — Ti < s — r < Si — Ti + e} 


we have C Wo © int(^) and Wi fl D = 0. 
(35) 


/c(s, r)h(s — r)h(r) dr 


>1^1 


h{s — 


It is not hard to show that 


T)h{T) dr 


- Y / l/(s-r)/(r)| dr 


for h G Lq{0,1) and (s,r) G Wi. Now we dehne sequences {h„}neN and {/n}neN 
by 

(= ) for ri < r < ri + £ 

hn{T~) '■= \ E ) for Si — ri < f < Si — ri + £ /n := / + r h^- 

\ 0 else 


Obviously ^ 0, fn^f in A, and due to the weak continuity of F (cf. Propo¬ 
sition]^ F{fn) F{f) in Y. With ( ]^ at hand it is easy to show that 

liminf ||F(hn)||Y > 0, 

n^oo 


which means that F{hn) 74 0 in A. Taking into account that the Frechet deriva¬ 
tive is compact and therefore completely continuous we have F'{f) —?■ 0 in Y. 
Together with formula (21), we then obtain 


F{fn) - F{f) = F{hn) + r F\f) K and thus F(/„) F{f). 
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The compactness of F\f) is an immediate consequence of the fact that (20) is a 
linear Fredholm integral operator with square integrable kernel. Such operators 
are always Hilbert-Schmidt operators and thus compact. The nonlinear operator 
F, however, is not compact, since the sequence {/njneN C X is bounded and 
weakly convergent to / but the associated sequence {F(/„)}„gN C Y cannot have 
a convergent subsequence. This completes the proof. □ 


Proof of Proposition |13[ Let {fn} C Xq be a sequence such that 



Passing to a subsequence if necessary, we may assume that ||F(/„)||y > e holds for 
all n. Due to the continuity of F(/), the pre-image of the closed set {||i/||y > e}, 


A:={feX : ||F(/)||y >4 

is a closed subset of X = 1). Now any accumulation point / of {/„} with 

respect to the weak topology in Xq (/ exists due to the stabilizing properties of 
TZ{f)) is a strong accumulation point of {fn} in X and hence contained in A^. 
Taking a subsequence, again denoted by {fn}, such that ^ / in Xq we thus 
have fn^ finX and F(/„) ^ F{f) in Y with ||F(/)||y > £. This yields 


Fjfn) - \F{fn)\F'l’' F{f) - \F{f)\e^^^ . 

Il^(/n)i|y ||X(/)||v 


and 


\fn\ I/I in X. 


In combination with the weak lower semicontinuity of F(/) in Xq we therefore 
obtain 


< limmfrf};!^{fn 

’ n^oo ’ 

and the proof is complete. 




□ 
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